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removes all dangerous terms which contain divergences. 
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§1. Introduction 

Recent experimental breakthrough of manufacturing Bose- Einstein condensates (BEC) 
in magnetically trapped dilute gases of alkali-metal atoms has created new opportunities to 
test many-body theories in well-controlled experimental conditions. Although the phenom- 
ena can be considered as a spectacular realization of an old prediction of quantum statistical 
mechanics as applied to ideal bose gases , it is known that mutual interactions of particles 
still play important roles to determine characteristic properties of the systems. 

In theoretical treatments the condensate is usually described by the dilute gas approxi- 
mation in which the pair-wise short range interaction is fully accommodated by the use of 
the effective interaction, usually referred to as the pseudo-potential, defined by 

V,^,{r) ^ ^^5{r) ^ g5{r) , (1-1) 

where a is the S'-wave scattering length of the pair interaction and m is the mass of the 
particle, while the many-body collective dynamics of the condensate is described by the 
Gross-Pitaevskii equation^), also known as the non-linear Schrodinger equation, for a clas- 
sical field ${r,t)] 



ih 



dt 
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rv 

2m 



+ gmr,t)\ 



^{r,t). (1-2) 



This equation can be derived from the Heisenberg equation of motion of the quantum field 
iP{r,t) 

in'-^ - [^(-'^)'^] - -'^^^ + 9i^Hr,mr,t) (1.3) 
by the replacement 

V'(r,i)^^(r,i) = (V'(r,i))eoh. 

where the expectation value is taken with a coherent state of the particle creation and 
annihilation operator of the lowest single particle state. This approximation is also known 
as the tree approximation in the language of the loop expansion^) and ignores quantum 
fluctuations around the mean fleld (?(r, t) . It is the purpose of the present work to investigate 
the effect of quantum corrections to the mean fleld by the variational method which we have 
developed for describing dynamics of quantum flelds. 

When the effective interaction of the form (1-1) is used to calculate the effect of the 
quantum fluctuations, a special care is needed to avoid (ultraviolet) divergences which arises 
due to the singular nature of the potential which causes coupling of all momentum modes. 
^) " It has been shown by Lee, Huang and Yang that such divergences can be eliminated 
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systematically by the proper use of the pseudo-potential defined by 



Vps.{r) = gS{r)—r = g5{r) + g5{r)r 



d_ 
dr 



(1-4) 



The difference, 



d_ 

dr ' 



(1-5) 



which has only vanishing contribution in the first order, generates divergent terms in higher 
order and it is precisely these divergent terms which cancel the divergences which are con- 



Physically this procedure is understood as the elimination of the double counting in the 
usual perturbation expansion in terms of the bare two particle interaction Vo(r): since the 
pair scattering (in vacuum) is already fully accounted for by the scattering length a, one 
should not include, when performing many body calculation, contributions from repeated 
use of the cfi^cctivc interaction (1-1) in the two particle scattering channel.*^ For example, 
the contributions from the diagram in Fig. 1 (a) should not be included in the many body 
calculation, because it is a part of the diagrams Fig. 1 (b) which is already included in the 
lowest order diagrams Fig. 1 (c) in terms of \4ff. • The problem of ultraviolet divergence can of 
course be avoided if one starts from a bare two-body interaction of finite range; this approach 
would however involve much more elaborate calculations even for a static system. 

Elimination of the ultraviolet divergences due to the singular potential looks technically 
similar to the renormalization of the effective local couplings in quantum field theories to 
eliminate ultraviolet divergences. ^^-^ This analogy was exploited further by Toyoda and 
more recently by Braaten and Nicto. However, here we choose to remove all divergences by 
simply avoiding double countings of pair scatterings, by inspection, based on the observation 
that the effective interaction (1-1) should be used only in the first Born "approximation" to 
describe two particle scattering. 

This docs not mean, however, that we should discard the quantum fluctuations entirely: 
the many-body interaction such as Fig. 2 should be still retained in the calculation. For 
example, the well-known Bogoliubov theory contains the effect of quantum fluctuations 
as zero-point oscillations of quasi-particle modes; the contribution to the energy density of 
this fluctuation would diverge if one uses the potential of the form (1-1). This divergence 
can be eliminated by the above procedure leaving behind non-analytic a^^^ correction which 
originates from the sum of genuine many-body interaction contributions. 

*^ Although the procedure to eliminate the ultraviolet divergence found in the literature '^'^^ is very 
suggestive, this interpretation of the subtraction of the double counting terms was never explicitly mentioned 
elsewhere to the best of our knowledge. 



tained in the higher order quantum corrections of V^s. • ^■^ 
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Fig. 1. The one- loop diagram depicted in (a) should not be included in the calculation since it is 
a part of (b) which is already included in the first order (tree) diagram (c). 




(a) (b) 
Fig. 2. Some diagrams which cannot be reduced to the pair interaction like Fig. 1 (b) and therefore 
should be included in the calculation. 

In this paper we present the basic formulation of the variational method and derive 
the equations of motion which generalizes the Gross-Pitaevskii mean field equation (1-2) 
to include the effect of quantum fluctuations in terms of the reduced density matrix. The 
resultant equations are similar to those known as the Hartree-Fock-Bogoliubov equations for 
fermions in nuclear many-boby theory and are identical, although slightly disguised, to 
the ones obtained for interacting bose-systems by Blaizot and Ripka^^^, and more recently 
by Griffin . The general static solutions of these equations are given in terms of the mode 
functions. We then perform an explicit calculation for a uniform condensate and show how 
the divergent terms appear by the use of the effective potential (1-1). By the inspection 
of the origin of these divergences in perturbation series we will demonstrate that all these 
divergences can be eliminated by the systematic removal of the terms which correspond to 
the double counting of pair scattering. By this method we obtain a gap equation free of 
divergence, (6-20), a new result obtained in this paper: To the best of our knowledge, this 
equation has never been written explicitly in the literature. In the forthcoming paper, we 
shall apply the present method to compute the density response function of the system. 

The rest of the paper is organized as follows. In the next section, we construct the 
Gaussian variational wave functional for condensates in a general trapping external potential. 
In section 3, we derive the equations of motion which consist of a generahzed Gross-Pitaevskii 
equation for the condensate and the Liouville equation for the reduced density matrix. Our 
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generalized Gross-Pitaevskii equations consist of two coupled non-linear equations: one is 
the equation of motion for <P{r,t) similar to the Gross-Pitaevskii equation (1-2), but it 
also contains the mean field terms associated with the fiuctuations, and the other is the 
equation of motion for the fiuctuations, written in the form of the Liouville equation for 
the reduced density matrix consisting of fiuctuations, these equations are often called the 
time-dependent Hartree-Bogoliubov equation. In section 4, we show that these equations 
can be easily extended to finite temperatures using the Gaussian density matrix which is 
determined self-consistently. In section 5, we express the static solutions of these equations 
in terms of the mode functions. In section 6 we examine the static solution in uniform 
system in detail and compare the results to the well-known Bogoliubov's mean field theory. 
We illustrate how the divergence can be eliminated by the removal of the double counting 
terms in the diagrams of perturbative series. A brief summary is given in section 7. 



§2. Gaussian wave functional 

In this and following sections we shall apply the variational method developed earlier 
in the functional Schrodinger picture to derive an extended form of the Gross-Pitaevskii 
equation which contain the effect of quantum fluctuations in a minimal fashion. 

We consider a non-relativistic gas of bose particles of mass m interacting via contact 
interaction in the form of (1-1). The dynamics of system is described by the following 
second quantized Hamiltonian: 

+^-i;\x)i;\x)i;{x)i;ix)^ . (2-1) 

Here, T4xt represents an external confining potential provided by the magnetic trap and the 
mutual interaction of particle is given by the effective interaction of the form (1-1) . The 
flelds ip{x) are quantized by the commutation relation : [ ip{x) , il^\y) ] = 5^{x — y) and 
the other combinations are equal to 0. Hereafter, we take h — 1. 

In order to apply the Gaussian wave functional method, which has been developed in 
our previous work^\ Schrodinger picture, we need to rewrite the Hamiltonian in terms of 
canonical coordinates and their associate momenta. There is some earlier attempt to deflne 
such coordinates but we choose here the representation which diagonalize the bilinear 
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part of the Hamiltonian. 



The field operator is decomposed as 



tp{x) =J2Va{x)aa 



(2-3) 



where ipa{x) is the normahzed eigen wave function of free bosons in external potential. We 
then introduce the following Hermitian operators : 

1 



Qa 



V'2'^a 

Pa^-^/§{a^-ai). (2-4) 
With / dx(p*^{x)ip()[x) = da,i3 these operators satisfy the canonical commutation relations: 

[Qa : ] = i^a,p , 

an analogue to that of the position and momentum operator in quantum mechanics, and 
convert the free Hamiltonian into a sum of uncoupled harmonic oscillators: 



(2-5) 



In the Q-representation, corresponding to the coordinate representation in quantum mechan- 
ics, the field momentum operator is expressed by the functional differential: Pa — —'^j^- 
The vacuum state of the Hamiltonian Hq is expressed as a product of Gaussian: 



-2 IZ^"<3a 



~ J d'xd'y<P^{x)Wo{x,y)<P{y] 



'Z^o[0] = (</'|O)=AAoexp 
= A/o exp 

where in the last line we have used the complex 0-field defined by* 

0(a;) =J2Va(x)Qa 

a 

and the kernel of the integral is defined by 

Wo{x,y) = J2v»{x)uJa^pUy) ■ 



(2-6) 



(2-7) 



*^ Do not confuse 4'{x) with the original bose field 'ip{x): it obeys the commutation relation, 
[(^(x), 7T^{y)] = i6{x — y) where the 7r-field is defined by 7r(a;) = (pa{x)Pa- 



6 



Indeed, this state will be annihilated when we operate the particle annihilation operator aa 
from the left for any a. In the uniform system (Kxt. = Vq) : (Pa{x) ~^ 
k'^/2m + Vo, and Wo{x,y) S{x - ?/)(-V|;/2m + I/q)- 

Our time- dependent variational wave functional in the functional Schrodinger picture is 
written as a straightforward generalization of the time-dependent Gaussian wave function in 
quantum mechanics: 

m = A^exp [i(7f 10 - 0) - (0 - + iE\<f> - 0)] (2-8) 

where we have used abbreviate notations for integrals such as 

(7f|0 - 4>) = E^aW {Qa - Qait)) = / te*(a;,i) - 4>{x,t)) , 

a,l3 ■' 

for K — G~^, U. Here two (complex) functions, 

a a 

and two kernels, 

G-\x,y-t) = Y.Mx)G:^^{t)^Uy) 

IJ{x,y;t) = Y,^a{x)UaAt)^Uy) 

are to be considered as time-dependent variational parameters of the wave functional. 

The expectation value of a composite operator tt) is given by the functional integral 

{o) = -^^^m 

which becomes a functional of the variational parameters (j){x,t), 7r{x,t), G{x,y,t) and 
U{x, y, t). For example, the expectation values of the field (plx) and its canonical momentum 
7r(a;) are given by 

{<f>{x)) ^ 4>{x,t) (2-9) 

{n{x)) ^ n{x,t) (2-10) 

respectively, and 

(V'(^)) = E J^Vaix) (Q^it) + -P„(t)) ^ ${X, t) , (2-11) 

ki'Kx)) = E (Oa W - ;f ^ r (a., . (2-12) 
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The nonvanishing expectation value of the particle creation and annihilation operators im- 
plies that our wave functional is not an eigen state of the particle number. Our Gaussian 
state corresponds to a squeezed state, or a generalized coherent state, in quantum optics. 
In order to describe field fluctuations it is convenient to introduce 

^{x,t) = ij{x)-^{x,t) . (2-13) 

For example, we obtain 

i^H^my)) = my, t) + {i^\x, t)^{y, t)) (2-i4) 

and so on, where {'ip\x,t)'ip{y,t)) is the two-point correlation function at time t and rep- 
resents the quantum fluctuation. The fluctuation is determined by the width parameter 
of the Gaussian wave functional. The explicit forms may be obtained from the following 
relations: 



= (2-15) 

1 

-( 

4 



{P^Pp) = iG-i(t) + 4 {SGE)^^^ it) (2-16) 



(4^/3) = -2 iGIJ)^^f, (t) + \i5^,f, (2-17) 
{P^Qp) = -2 {EG)^^^ it) - (2-18) 

where = Qa ~ Qa and Pa = Pa — Pa- For example, 

{4){x)4)^ (y)) = J2^»i^){QocQi3)^}i{y) 

^J2^aix)Ga,p{t)(p*p{y) = G{x,y;t) 

and similar relations can be obtained for {(j)\x)Ti{y)) etc. However, we do not need such 
explicit forms in the following discussion. 

With the Gaussian form of the variational wave functional, expectation value of all higher 
products of the field operators are expressed in terms of the center (the mean field) and the 
width (the two-point correlation function) of the Gaussian. This feature, characteristic of 
the Gaussian Ansatz of the fluctuations, enables one to truncate otherwise infinite set of 
coupled equations of the motion for the expectation values of the products of the quantum 
fields into a set of coupled equations of motion for the mean field ^{x, t) and the two-point 
functions of fiuctuations G{x,y;t). 
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^3. Equations of motion 



The equations of motion can be derived from the time-dependent variational principle as 
applied to the action: 

S = I dmi--H + piN\^) , (3-1) 

where /i is the chemical potential of the bosons introduced to enforce that the total number 
of the bosons N — (WlN'l'I^), defined with N = J dxilj\x)ilj{x), is conserved. Minimization 
of (3-1) with respect to the variational parameters leads to the equations of motion of (j){x, t) 
and 7t{x, t), G~^{x, y; t) and S{x, y; t) which can be rewritten for the mean field and the 
fluctuations. 

Alternatively, one can derive the equations of motion written in terms of the mean field 
and fiuctuation directly from: = 6{H')/S{^*) and i^* = -5{H')/5<P. Here, {H') is the 
expectation value of the effective Hamiltonian H' — H — jiN . Then one finds: 

iS{x,t) = -^V^<P{x,t) + (KxtN - ^J)${x,t) + g\<P{x,t)\^<P{x,t) 

+2g{'^'^{x)i^{x))${x, t) + g{i^{x)'^{x))$*{x, t) , (3-2) 

which is a generalized form of the Gross-Pitaevskii equation (1-2) including the effect of 
fluctuations. 

Our next task is to derive equations of motion for the two-point functions: {■ijj'^{x)'i^{y)) 
and {il'{x)i/j{y)). This can be done from the equations of motion of the kernel functions 
G{x,y,t) and E{x,y,t). For this purpose we introduce the 2x2 reduced density matrix 
deflned as 

M{x,y,t) + -S%x-y)^[ . . /iu\i(\\]- ^^'^^ 

Then one can show that the equation of motion for the reduced density matrix Ai is expressed 
in the form of the Liouville equation: 

iM^[n , M] , (3-4) 

where the Hamiltonian density 7i is defined by 



with 



n(x,y,t) = S''(x-y) { \\ \\\ (3-5) 



W{x) = -^V^ + V^x) - // + 2g\^{x, t)]' + 2g{i^\x)i^{x)) 

= W*(x) (3-6) 
A{x) = g{<P{x, tf + {^{x)^{x))) . (3-7) 
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As we shall show, this equation reproduces the quasiparticle excitations in the Bogoliubov 
theory when the term induced by the fluctuations are discarded. 

Equations (3-2) and (3-4) with (2-11), (3-3), (3-5), (3-6), and (3-7), form a closed set of 
equations which describe non-linear time evolution of the Bose-Einstein condensate and its 
fluctuations described by the Gaussian state. They are an extension of the Gross-Pitaevskii 
equation to include quantum fluctuations. 



\4. Inclusion of thermal fluctuations 



Foregoing derivation based on the pure state (2-8) can be extended to include thermal 
fluctuations by introducing statistical average with a Gaussian density matrix for mixed 
states. If the system is near equilibrium, it is equivalent to introduce the density operator 
in the form: 

D{t) = le-^^(*) , with Z = Tre-^^(*) , (4-1) 

Zj 

where ji = is the inverse temperature and H{t) is the (time-dependent) mean fleld 

Hamiltonian given by 



H{t) = I d?x{ -—^l^\x)V^il^{x) + (V;.t(x) - ii)ilj\x)tlj{x) 



{ip\x)ip\x))ip{x)ip{x) + {ip{x)ip{x))^\x)^\x) 



+4:{ij^{x)^{x))^\x)ij{x] 



(4-2) 



Here, {^j^ {x)ij{x)) = $*{x,t)^{x,t) + {^j^ {x)ij{x)) , {'^\x)ij^ (x)) = {x , t)$* {x , t) + 
{'ij)\x)-ilj\x)) and {tp{x)ip{x)) = ^{x,t)${x,t) + {ip{x)ip{x)). The mean field and the fiuc- 
tuations are now defined by 



^ix,t) = iijix)) = Tr(D(t)^(a;)) 
{^{x)^{y)) = Tr (D{t)^{x)^{y)) etc. 



(4-3) 
(4-4) 

One can derive the equation of motion from the Liouville equation of the non-equilibrium 
density matrix D, 



dD 



H,D 



For example, the equation of motion of the mean field can be obtained from 



dt 



iTi (^b{t)ijj{x)) = Tr ([//, d\ ij{x)) = Tr [^(a;), H 



(4-5) 



(4-6) 
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Upon insertion of the (4-2) to compute the commutator H,iIj{x) we reproduce the gener- 
ahzed Gross-Pitaevskii equation in the form as (3-2) which we have obtained for the pure 
state evolution. Similarly, by replacing ip^x) by '!jj{x)'ijj{y), etc. one can derive the equa- 
tion of motion for the reduced density matrix Ai defined by (3-3) with (4-4) which becomes 
precisely the same as the equation (3-4). 

We can therefore reinterpret that the equations, (3-2) and (3-4), also describe the non- 
equilibrium time evolution of the perturbation introduced in the system which has been 
in equilibrium at temperature T and the chemical potential /x. 

§5. Static solutions 

We first consider the case in which the mean field ^ and the fluctuation Ai do not depend 
on time, that is, the system is in the thermal equilibrium. In this case, the equation of motion 
for M. becomes 

[n , M ] = . (5-1) 

This implies that the Hamiltonian density H and the reduced density matrix M. can be 
simultaneously diagonalized. This is achieved explicitly by introducing the mode functions. 
The mode functions, u{x,t) and v{x,t), for the Hamiltonian matrix Ti are deflned by: 

(uix,t)\ ^ / Wix) A{x) \ (u{x,t)\ 
''\v{x,t)) \-A*{x) -W*{x))\v{x,t)) - ^ ^ 

For the stationary state, we write 

u{x, t) = e-'^^u{x) , v{x, t) = e-'^^v{x) . (5-3) 

The eigenvalue equation is then obtained as 

f W{x) A{x) \U^{x)\^^ fu^{x)\ 
\-A*{x) -W*{x)j\v^{x)) "UaN/ ' 
where a represents quantum numbers. By taking a complex conjugate of the above eigenvalue 
equation, it can be shown that the following equation is also satisfied : 
f W{x) A{x) \fvl{x)\^ fvl{x)\ 
\-A*{x) -W*{x)) \ul{x)) n<Wy ■ ^ ^ 

Namely, there exist negative energy solutions. The eigenvector (m, v) is also the eigenvector 
for the reduced density matrix M. because H, and M. commute : 

jd'y{M{x,y) + \5\x-y)) ([^j^j + ^2) ) • (5-6) 

These equations describe the coherent, adiabatic evolution of the system, while kinetic aspects (e.g. 
the growth of the condensate ' ) are beyond the scope of the present method. 
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where, by construction (3-3), the eigenvalues of M. are related to the quasi-particle distri- 
bution Ua by^'''^^ 

1 1 

fa = na + - , na= ^^^^ _ ^ ■ (5-7) 

Taking the complex conjugate of the eigenvalue equation for + 1/2, we also obtain 
Here the ortho-normalization relation is given by 

j dPx {ul{x)Ufi{x) - V*^{x)vp{x)) = ±6af3 , (5-9) 

where a > (< 0) and /3 > (< 0) corresponds to -|- (— ) on the right-hand side and it is 
implied that a > (< 0) represents the positive (negative) energy solution. The reduced 
density matrix can be expressed in terms of the mode functions u and v as 

M{x, y) + U\x - 2/) = E ( f """Ji ] Ua{x) + 1/2) «(2,), -<(?/) ) 

-([]:|^|)(/aN-l/2)(-^a(2/), «a(y))| . (5-10) 

Here, summation is taken only over the positive energy eigenstates. This result implies, 
from the definition (3-3), the matrix elements of the reduced density matrix are expressed 
in terms of the mode functions and the bose distribution function n„ : 



{i)\x)i){y)) = Y,[{naix) + l)vaix)vl{y) + na{x)ul{x)ua{y)] 
{tP{x)tp{y)) = + l)uaix)vl{y) + na{x)v*^{x)ua{y)] 

{'^\x)'^\y)) = Y.[Mx) + l)va{x)ul{y) + n^{x)ul{x)va{y)\ 
{ij{x)ij\y)) = ^[KN + l)«a(a;)<(y) + n^{x)vl{x)vM ■ (5-11) 

a>0 

The equations (3-2), (5-4) and (5-11) with (5-7) form a basic set of equations for the static 
Bose-Einstein condensates which need to be solved self- consistently. 

§6. Uniform condensate 

In this section, we consider a simple uniform system, neglecting external confining po- 
tential, V^xt = 0, in order to illustrate how our method works. In this case, the condensate 
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is spatially uniform: (I>{x) = and all mode functions are reduced to simple plane waves. 
The indices a are regarded as momenta k and the mode functions u and v may be written 
as 

Ua — ' Mk(x) = -^e'^-^Uk , Va — > Wk(x) = -^e'^'^Vi, , (6-1) 

where V is the volume of the system. Further, we take the continuum limit {V — > oo) with 
the prescription J2k — V/{27r)^ • J o?^k and Sap — (^(k — k'). 
The generahzed Gross-Pitaevskii equation now takes the form, 

(-/i + g\$o\' + 25(^^^))^o + 9{^^)^*o = , (6-2) 

where 

(V'V) = J ^[(rik + l)|^;k|' + nk|«kr] , (6-3) 
r d^k 

{^^) = J j^{2n^ + l)uy,v*^ (6-4) 
and the mode functions and their eigenvalues are easily obtained from (5-4) as 



l^k + ^k W k-^k ^* ,^ 



with 



= ^W^ - \A\^ , (6-6) 

where 

W^k = 1^ - + 2^(|^or + (V^^^)) , (6-7) 

^ = ^(^o' + (^^)) • (6-8) 

Inserting (6-4) into (6-8) and then using (6-5), we obtain the "gap equation" for the 
anomalous density: 

- = {^^) = ^0 + {U) 
9 



On the other hand, the particle number density is given by 

2 f d^k ( W^k-^k 
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where we have used (6-5) again. Equations, (6-2) and (6-9), together with the constraint 
(6-10) and with the quasi-particle distribution function (5-7), determine Z\, <fo and /i, self- 
consistently for a given value of particle density n and temperature T. 

We can show that the well-known Bogoliubov spectrum of the quasi-particle can be 
reproduced from (6-6) by removing all fluctuations. For simplicity, we consider the ground 
state, nk = for all fc 7^ and assume that is real. Then, by setting {ip'^'ip) = {tp''P) ~ 
{tp'^'ij)'^) = 0, we find 

l^ = \A\=g\<PoW (6-11) 

and 




where we have used 



E^. = VW^k - W = A 7^ + 2^l^oP , (6-12) 



= ^ - ^ + 2g\<P,\^ = 9\^o\' , (6-13) 

2m Zm 



which precisely coincides with the Bogoliubov dispersion relation for quasi-particle excita- 
tions. For small with 5* > 0, the phonon dispersion relation uj = ck is obtained, where 
the sound velocity c is 



(6-14) 
V m 

This is the well-known result of the Bogoliubov theory. 

The appearance of the phonon dispersion relation is a consequence of spontaneous break- 
down of the continuous U{1) symmetry of the Hamiltonian associated with the phase change 
of the bose field: '4>{x) e^^il>{x). The phonon can be regarded as a Nambu-Goldstone mode 
of the broken U{1) (or 0(2)) symmetry. 

Our formula (6-6) for the spectrum of the mode functions with the Gaussian fiuctuation, 
however, does not reproduce phonon spectrum, but instead exhibits a gap at k — > 0, violating 
Goldstone's theorem. This shortcoming is rather inherent in the Gaussian wave functional 
approach, but this defect can be remedied by computing the excitation spectrum. 

We now turn our discussion on the removal of divergences. The expression (6-4) for the 
fiuctuations in the reduced density matrix actually contain a divergent integral, and so does 
the gap equation (6-9). To see this expficitly, we observe that at large |k|. 
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(6-15) 



*^ According to the classification of Hohenberg and Martin our approximation corresponds to the 
conserving (or (^-derivable ) approximation, as opposed to the gapless approximation^^). The origin of the 
violation of the Goldstone theorem is traced back to the "Fock" term which may be suppressed by keeping 
the leading term in the 1/N expansion^") for a system with 0(N) symmetry. 
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Fig. 3. The divergent contribution of the off-diagonal component of the density matrix, 
so that 



Inserting this expression in (64), we find that the off-diagonal component of the density 
matrix becomes 

j-^(2n, + 1K< -/ (^(2». + 1)^ , (6-17) 
which indeed contains a hnearly divergent integral, 

This divergent integral is associated with the diagram shown in Fig. 3. It originates 
from a repeated use of the effective potential as we have mentioned in the introduction. 
This divergence should be removed in order to avoid the double counting of the perturbation 
series in terms of the bare potential. Therefore we define the renormalized density matrix 
by subtracting the divergent term: 

(^^)ren. = (^^) " (V'^)div. (6-19) 

This procedure is equivalent to replace the gap equation by: 

With this replacement and all other relations untouched, all divergences associated with 
double counting of the pair scatterings are removed from the previous equations. The "renor- 
malized gap equation" (6-20) is one of our new findings in this paper. 



§7. Summary 

In this paper we have formulated the variational method to describe the dynamics of 
the Bose-Einstein condentate interacting via short-range effective interaction. With the 
Gaussian variational wave functional we have derived the equations of motion which consist 
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of a generalized Gross-Pitaevskii equation for the condensate and a Liouville-von Neumann 
equation for the density matrix of the quantum and thermal fluctuations. The static solutions 
of these equations are examined in terms of the mode functions for the fluctuations. In the 
case of uniform condensate, these mode functions are reduced to simple plane waves, and 
we are able to perform the calculation explicitly. We have shown that the spectrum of the 
mode functions reduces to the Bogoliubov phonon spectrum when we neglect the effect of 
quantum fluctuations; it does not contain however phonon mode in the presence of quantum 
fluctuation. In these calcuations we have shown that the divergent integrals which arise due 
to the use of singular eflective interaction can be eliminated systematically by removing the 
repeated pair scattering terms to avoid the double counting of the diagrams in terms of the 
bare interaction. In the forthcoming paper, we will show how to apply the present method 
to compute the linear response of the system to external perturbation and to extract the 
phonon dispersion relation of the collective excitations in the sytem. 
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